#include "HY1C_out.h"
#include <stdlib.h>
#include <math.h>
#include "l12_proto.h"
#include "chl.h"


float chl_oc2(l2str *l2rec, float Rrs[])
{
    static int32_t  *w  = NULL;
    static float *a  = NULL;
    static int   ib1 = -1;
    static int   ib2 = -1;

    float rat;
    float Rrs1, Rrs2;
    float chl = chlbad;

    if (w == NULL) {
        w = l2rec->input->chloc2w;
        a = l2rec->input->chloc2c;
        if (w[0] < 0 || w[1] < 0) {
            HY1C_out("chl_oc2: algorithm coefficients not provided for this sensor.\n");
            exit(1);
	}
        ib1 = bindex_get(w[0]);
        ib2 = bindex_get(w[1]);
        if (ib1 < 0 || ib2 < 0) {
            HY1C_out("chl_oc2: incompatible sensor wavelengths for this algorithm\n");
            exit(1);
        }
    }

    Rrs1 = Rrs[ib1];
    Rrs2 = Rrs[ib2];

    if (Rrs1 > 0.0 && Rrs2 > 0.0 && Rrs1/Rrs2 <= 10.0) {
        rat = Rrs1/Rrs2;
        if (rat > minrat && rat < maxrat) { 
            rat = log10(rat);
            chl = (float) 
               pow(10.0,(a[0]+rat*(a[1]+rat*(a[2]+rat*(a[3]+rat*a[4])))));
            chl = (chl > chlmin ? chl : chlmin);
            chl = (chl < chlmax ? chl : chlmax);
	}
    }

    return (chl);
}    


float chl_oc3(l2str *l2rec, float Rrs[])
{
    static int32_t  *w  = NULL;
    static float *a  = NULL;
    static int   ib1 = -1;
    static int   ib2 = -1;
    static int   ib3 = -1;

    float rat, minRrs;
    float Rrs1, Rrs2, Rrs3;
    float chl = chlbad;

    if (w == NULL) {
        w = l2rec->input->chloc3w;
        a = l2rec->input->chloc3c;
        if (w[0] < 0 || w[1] < 0 || w[2] < 0) {
            HY1C_out("chl_oc3: algorithm coefficients not provided for this sensor.\n");
            exit(1);
	}
        ib1 = bindex_get(w[0]);/*
        ib2 = bindex_get(w[1]);
        ib3 = bindex_get(w[2]);*/
        ib1 =1;//443
        ib2 = 2;//490
        ib3 = 4;//565
        if (ib1 < 0 || ib2 < 0 || ib3 < 0) {
            HY1C_out("chl_oc3: incompatible sensor wavelengths for this algorithm\n");
            exit(1);
        }
    }

    Rrs1 = Rrs[ib1];
    Rrs2 = Rrs[ib2];
    Rrs3 = Rrs[ib3];

    minRrs = MIN(Rrs1,Rrs2);

    if (Rrs3 > 0.0 && Rrs2 > 0.0  && minRrs > -0.001) {
        rat = MAX(Rrs1,Rrs2)/Rrs3;
        if (rat > minrat && rat < maxrat) { 
            rat = log10(rat);
            chl = (float) 
               pow(10.0,(a[0]+rat*(a[1]+rat*(a[2]+rat*(a[3]+rat*a[4])))));
            chl = (chl > chlmin ? chl : chlmin);
            chl = (chl < chlmax ? chl : chlmax);
	}
    }

    return (chl);
}    


float chl_oc3c(l2str *l2rec, float Rrs[])
{
    static float a[]  = {0.2515,-2.3798,1.5823,-0.6372,-0.5692};
    static int   ib1 = -1;
    static int   ib2 = -1;
    static int   ib3 = -1;

    float rat, minRrs;
    float Rrs1, Rrs2, Rrs3;
    float chl = chlbad;

    if (ib1 < 0) {
        ib1 = bindex_get(443);
        ib2 = bindex_get(490);
        ib3 = bindex_get(545);
        if (ib3 < 0) ib3 = bindex_get(550);
        if (ib3 < 0) ib3 = bindex_get(555);
        if (ib3 < 0) ib3 = bindex_get(560);

        if (ib1 < 0 || ib2 < 0 || ib3 < 0) {
            HY1C_out("chl_oc3: incompatible sensor wavelengths for this algorithm\n");
            exit(1);
        }
    }

    Rrs1 = Rrs[ib1];
    Rrs2 = Rrs[ib2];
    Rrs3 = Rrs[ib3];

    minRrs = MIN(Rrs1,Rrs2);

    if (Rrs3 > 0.0 && Rrs2 > 0.0  && minRrs > -0.001) {
        Rrs3 = conv_rrs_to_555(Rrs3,l2rec->fwave[ib3]);
        rat = MAX(Rrs1,Rrs2)/Rrs3;
        if (rat > minrat && rat < maxrat) { 
            rat = log10(rat);
            chl = (float) 
               pow(10.0,(a[0]+rat*(a[1]+rat*(a[2]+rat*(a[3]+rat*a[4])))));
            chl = (chl > chlmin ? chl : chlmin);
            chl = (chl < chlmax ? chl : chlmax);
	}
    }

    return (chl);
}    


float chl_oc4(l2str *l2rec, float Rrs[])
{
    static int32_t  *w  = NULL;
    static float *a  = NULL;
    static int   ib1 = -1;
    static int   ib2 = -1;
    static int   ib3 = -1;
    static int   ib4 = -1;

    float rat, minRrs;
    float Rrs1, Rrs2, Rrs3, Rrs4;
    float chl = chlbad;

    if (w == NULL) {
        w = l2rec->input->chloc4w;
        a = l2rec->input->chloc4c;
        if (w[0] < 0 || w[1] < 0 || w[2] < 0 || w[3] < 0) {
            HY1C_out("chl_oc4: algorithm coefficients not provided for this sensor.\n");
            exit(1);
	}
        ib1 = bindex_get(w[0]);
        ib2 = bindex_get(w[1]);
        ib3 = bindex_get(w[2]);
        ib4 = bindex_get(w[3]);
        if (ib1 < 0 || ib2 < 0 || ib3 < 0 || ib4 < 0) {
            HY1C_out("chl_oc4: incompatible sensor wavelengths for this algorithm\n");
            exit(1);
        }
    }

    Rrs1 = Rrs[ib1];
    Rrs2 = Rrs[ib2];
    Rrs3 = Rrs[ib3];
    Rrs4 = Rrs[ib4];

    minRrs = MIN(Rrs1,Rrs2);

    if (Rrs4 > 0.0 && Rrs3 > 0.0 && (Rrs2 > 0.0 || Rrs1*Rrs2 > 0.0) && minRrs > -0.001) {
        rat = MAX(MAX(Rrs1,Rrs2),Rrs3)/Rrs4;
        if (rat > minrat && rat < maxrat) { 
            rat = log10(rat);
            chl = (float) 
               pow(10.0,(a[0]+rat*(a[1]+rat*(a[2]+rat*(a[3]+rat*a[4])))));
            chl = (chl > chlmin ? chl : chlmin);
            chl = (chl < chlmax ? chl : chlmax);
	}
    }

    return (chl);
}    

 
float chl_hu(l2str *l2rec, float Rrs[])
{
    static float  w[]  = {443,555,670};
    static float  a[]  = {-0.4909,191.6590};
    static int   ib1 = -1;
    static int   ib2 = -1;
    static int   ib3 = -1;

    float base;
    float Rrs1, Rrs2, Rrs3;
    float chl  = chlbad;

    if (ib1 == -1) {
        ib1 = bindex_get(443);
        ib2 = bindex_get(545);
        if (ib2 < 0) ib2 = bindex_get(550);
        if (ib2 < 0) ib2 = bindex_get(555);
        if (ib2 < 0) ib2 = bindex_get(560);
        ib3 = bindex_get(670);
        if (ib1 < 0 || ib2 < 0 || ib3 < 0) {
            HY1C_out("chl_hu: incompatible sensor wavelengths for this algorithm\n");
            exit(1);
        }
        // reset to actual center wavelengths
	//w[0] = l2rec->fwave[ib1];
        //w[1] = l2rec->fwave[ib2];
        //w[2] = l2rec->fwave[ib3];
    }

    Rrs1 = Rrs[ib1];
    Rrs2 = Rrs[ib2];
    Rrs3 = Rrs[ib3];

    if (Rrs3 > BAD_FLT+1 && Rrs2 > 0.0 && Rrs1 > 0.0) {
        Rrs2 = conv_rrs_to_555(Rrs2,l2rec->fwave[ib2]);
        base = MIN(Rrs2 + (w[1]-w[0])/(w[2]-w[0])*(Rrs3 - Rrs1),0.0);
        if (base <= 0.0) { 
            chl = (float) pow(10.0,a[0] + a[1]*base);
            chl = (chl > chlmin ? chl : chlmin);
            chl = (chl < chlmax ? chl : chlmax);
	}
    }

    return (chl);
}    

 
float chl_oci(l2str *l2rec, float Rrs[])
{
    static float t1 = 0.25;
    static float t2 = 0.30;

    float chl1 = chlbad;
    float chl2 = chlbad;
    float chl  = chlbad;

    chl1 = chl_hu(l2rec,Rrs);
    if (chl1 <= t1)
        chl = chl1;
    else {
        chl2 = get_default_chl(l2rec,Rrs);
        if (chl2 > 0.0) {
            if (chl1 >= t2)  
                chl = chl2;
            else {
	        chl = chl1 * (t2-chl1)/(t2-t1)
		    + chl2 * (chl1-t1)/(t2-t1);
                //chl = MIN(MAX(chl,t2),t1);
            }
	}
    }

    return (chl);
}    

 
float chl_clark(l2str *l2rec,float nLw[])
{
    static int32_t  *w  = NULL;
    static float *a  = NULL;
    static int   ib1 = -1;
    static int   ib2 = -1;
    static int   ib3 = -1;

    float rat;
    float nLw1, nLw2, nLw3;
    float chl = chlbad;

    if (w == NULL) {
        if (l2rec->input->outband_opt >= 2) {
            HY1C_out("chl_clark: this algorithm requires full-band nLw\n");
            HY1C_out("chl_clark: you must disable the out-of-band correction for nLw\n");
            HY1C_out("chl_clark: set outband_opt = 0\n");
            exit(1);
	}
        w = l2rec->input->chlclarkw;
        a = l2rec->input->chlclarkc;
        if (w[0] < 0 || w[1] < 0 || w[2] < 0) {
            HY1C_out("chl_clark: algorithm coefficients not provided for this sensor.\n");
            exit(1);
	}
        ib1 = bindex_get(w[0]);
        ib2 = bindex_get(w[1]);
        ib3 = bindex_get(w[2]);
        if (ib1 < 0 || ib2 < 0 || ib3 < 0) {
            HY1C_out("chl_clark: incompatible sensor wavelengths for this algorithm\n");
            exit(1);
        }
    }

    nLw1 = nLw[ib1];
    nLw2 = nLw[ib2];
    nLw3 = nLw[ib3];

    if (nLw1 > 0.0 && nLw2 > 0.0 && nLw3 > 0.0) {
        if (ib1 == ib2)
            rat = log10(nLw1/nLw3);
        else 
            rat = log10((nLw1+nLw2)/nLw3);
        chl = (float) 
            pow(10.0,(a[0] + rat*(a[1] + rat*(a[2] + rat*(a[3] + rat*(a[4] + rat*a[5]))))));
        chl = (chl > chlmin ? chl : chlmin);
        chl = (chl < chlmax ? chl : chlmax);
    }

    return(chl);
}    


float get_default_chl(l2str *l2rec, float Rrs[])
{
    static int   firstCall = 1;
    float chl;

    if (firstCall) {
        firstCall = 0;
    }

    chl = chlbad;

    switch (l2rec->sensorID) {
	case SEAWIFS:
	case OCTS:
	case OCM1:
	case OCM2:
	case MOS:
	case MERIS:
	case MERISL2:
            chl = chl_oc4(l2rec,Rrs);
            break;
	case MODIST:
	case MODISA:
	case HMODIST:
	case HMODISA:
	case CZCS:
	case OSMI:
        case VIIRS:
        case OCRVC:
            chl = chl_oc3(l2rec,Rrs);
            break;
	default:
        chl = chl_oc3(l2rec,Rrs);
        //20180824
//	    HY1C_out("%s Line %d: need a default chlorophyll algorithm for this sensor\n",
//                __FILE__,__LINE__);
 //           exit(1);
            break;
    }

    return(chl);
}


void get_chl(l2str *l2rec, int prodnum, float prod[])
{
    int32_t  ip;
    int32_t  ipb;

    static int   firstCall = 1;

    if (firstCall) {
        firstCall = 0;
    }

    /*                                                      */
    /* Compute desired products at each pixel               */
    /*                                                      */
    for (ip=0; ip<l2rec->npix; ip++) {

        ipb = ip*NBANDS;

        switch (prodnum) {

	    case DEFAULT_CHL:
                prod[ip] = get_default_chl(l2rec,&l2rec->Rrs[ipb]);
                break;
	    case CAT_logchl:
                prod[ip] = get_default_chl(l2rec,&l2rec->Rrs[ipb]);
                if (prod[ip] > 0.0) 
                    prod[ip] = log(prod[ip]);
                else
                    prod[ip] = BAD_FLT;
                break;
	    case CAT_log10chl:
                prod[ip] = get_default_chl(l2rec,&l2rec->Rrs[ipb]);
                if (prod[ip] > 0.0) 
                    prod[ip] = log10(prod[ip]);
                else
                    prod[ip] = BAD_FLT;
                break;
  	    case CAT_chl_oc2:
                prod[ip] = chl_oc2(l2rec,&l2rec->Rrs[ipb]);
                break;
  	    case CAT_chl_oc3:
                prod[ip] = chl_oc3(l2rec,&l2rec->Rrs[ipb]);
                break;
  	    case CAT_chl_oc3c:
                prod[ip] = chl_oc3c(l2rec,&l2rec->Rrs[ipb]);
                break;
  	    case CAT_chl_oc4:
                prod[ip] = chl_oc4(l2rec,&l2rec->Rrs[ipb]);
                break;
  	    case CAT_chl_hu:
                prod[ip] = chl_hu(l2rec,&l2rec->Rrs[ipb]);
                break;
  	    case CAT_chl_oci:
                prod[ip] = chl_oci(l2rec,&l2rec->Rrs[ipb]);
                break;
  	    case CAT_chl_clark:
                prod[ip] = chl_clark(l2rec,&l2rec->nLw[ipb]);
                break;
            default:
                HY1C_out("Error: %s : Unknown product specifier: %d\n",__FILE__,prodnum);
                exit(FATAL_ERROR);
                break;
        }

        if (prod[ip] == chlbad) 
            l2rec->flags[ip] |= PRODFAIL;
    }
}

